A Counterfactual Simulation Model of Causal Judgment
Tobias Gerstenberg
June 07, 2020
1 Setup
1.1 Load packages
library("knitr") # for knitting stuff
library("kableExtra") # for markdown tables
library("lme4") # for linear mixed effects models
library("broom") # for tidying model results
library("brms") # Bayesian regression with Stan
library("corrr") # for tidy correlation matrix
library("xtable") # for latex tables
library("jsonlite") # for reading json files
library("janitor") # cleans stuff
library("Hmisc") # bootstrapped confidence intervals
library("tidybayes") # for Bayesian data analysis
library("png") # adding pngs to images
library("grid") # functions for dealing with images
library("plotly") # 3D scatter plot
library("egg") # for geom_custom
# library("ggtern") # for ternary plots
library("patchwork") # for figure panels
library("tidyverse") # for wrangling, plotting, etc. 1.2 Helper functions
# setting some knitr options
opts_chunk$set(comment = "",
fig.show = "hold")
# set the ggplot theme
theme_set(theme_classic())
# function for printing out html or latex tables
print_table = function(data, format = "html", digits = 2){
if(format == "html"){
data %>%
kable(digits = digits) %>%
kable_styling()
}else if(format == "latex"){
data %>%
xtable(digits = digits) %>%
print(include.rownames = F,
booktabs = T)
}
}
# root mean squared error
rmse = function(x, y){
return(sqrt(mean((x - y)^2)))
}
actual_counterfactual_plot = function(x, ylabel) {
p = ggplot(data = x,
mapping = aes(x = clipindex,
y = value,
fill = colorindex,
group = model)) +
stat_summary(geom = "bar",
fun = "mean",
color = "black",
position = position_dodge(0.7),
width = 0.7) +
stat_summary(geom = "linerange",
fun.data = "mean_cl_boot",
position = position_dodge(0.7),
size = 1) +
geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
shape = 21,
color = "black",
size = 1,
alpha = 0.2) +
scale_fill_manual(values = c("red2", "green2", "black")) +
facet_grid(index_actual ~ index_counterfactual) +
geom_text(data = df.text %>%
drop_na(label),
mapping = aes(x = x, y = y, label = as.character(label)),
size = 6,
position = position_dodge(width = .7),
hjust = 0.5) +
coord_cartesian(xlim = c(0.5, 2.5),
ylim = c(-5, 105),
expand = F,
clip = "off") +
labs(y = ylabel) +
theme(text = element_text(size = 20),
legend.position = "none",
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
panel.spacing.y = unit(0.8, "cm"),
plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
axis.title.x = element_blank(),
panel.grid = element_blank(),
strip.background = element_blank(),
panel.border = element_rect(colour = "black", fill = NA))
# print(p)
}
actual_counterfactual_threeballs_plot = function(x, ylabel) {
p = ggplot(data = x,
mapping = aes(x = ball,
y = value,
group = model,
fill = colorindex)) +
stat_summary(fun = "mean",
geom = "bar",
color = "black",
position = position_dodge(0.9),
width = 0.9) +
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
color = "black",
position = position_dodge(0.9)) +
geom_point(data = x %>%
mutate(value = ifelse(model == "rating", value, NA)),
position = position_jitterdodge(jitter.width = 0.3,
jitter.height = 0,
dodge.width = 0.9),
shape = 21,
color = "black",
size = 1,
alpha = 0.15) +
facet_wrap(~clip, ncol = 8) +
geom_text(data = df.text,
mapping = aes(x = x,
y = y,
label = as.character(label)),
size = 5,
position = position_dodge(width = .9),
hjust = 0.5,
fontface = "bold") +
coord_cartesian(xlim = c(0.5, 2.5),
ylim = c(-5, 120),
expand = F,
clip = "off") +
scale_y_continuous(breaks = seq(0, 100, 25)) +
scale_fill_manual(values = c("red2", "green2", "black")) +
scale_color_manual(values = c("red2", "green2", "black")) +
labs(y = ylabel) +
theme(text = element_text(size = 16),
legend.position = "none",
axis.title.x = element_blank(),
panel.spacing.y = unit(0.3, "cm"),
panel.grid = element_blank(),
strip.background = element_blank(),
strip.text = element_blank(),
panel.border = element_rect(colour = "black", fill = NA))
}2 Experiment 1: 2 balls
2.1 Clips
2.2 Read in data
# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv",
header = T,
stringsAsFactors = F)
# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
header = T,
stringsAsFactors = F)
# Information about each clip
df.exp1.clipinfo = tibble(clip = 1:18,
outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 1),
outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1,
0, 1, 1, 0, 0, 1, 0, 1, 1),
index_actual = rep(c("actual miss",
"actual close",
"actual hit"), each = 6),
index_counterfactual = rep(rep(c("counterfactual miss",
"counterfactual close",
"counterfactual hit"),
each = 2),
3)) %>%
mutate(index_actual = factor(index_actual, levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
# Structure data
df.exp1.causal.long = df.exp1.causal$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:18, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
arrange(participant, clip, rating) %>%
left_join(df.exp1.clipinfo,
by = "clip")
df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:18, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
arrange(participant, clip, rating) %>%
left_join(df.exp1.clipinfo,
by = "clip")2.2.1 Demographics
# counterfactual condition
df.exp1.counterfactual %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time"),
sep = ",") %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
select(gender:time) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 41 | 32 | 8.8 | 19 | 7.5 | 3.12 |
# causal condition
df.exp1.causal %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time"),
sep = ",") %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
select(gender:time) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 41 | 33 | 10.3 | 21 | 9.65 | 6.58 |
2.3 Read in model predictions
2.3.1 Counterfactual condition
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "2ball")]
df.exp1.model = files %>%
map_dfr(~ read.csv(str_c(path, .))) %>%
set_names(c("clip", "noise", "prediction")) %>%
left_join(df.exp1.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))
# calculate mean counterfactual judgments
df.exp1.counterfactual.means = df.exp1.counterfactual.long %>%
group_by(clip) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup() %>%
left_join(df.exp1.clipinfo, by = "clip")`summarise()` ungrouping output (override with `.groups` argument)# find noisy simulation model that best predicts the mean counterfactual judgments by
# minimizing the sum of squared errors
df.exp1.counterfactual.model = df.exp1.model %>%
group_by(noise) %>%
nest() %>%
mutate(sse = map(data,
~ sum((.x$prediction*100 -
df.exp1.counterfactual.means$rating_mean) ^ 2)),
sse = unlist(sse)) %>%
ungroup() %>%
filter(sse == min(sse)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)2.3.2 Causal condition
df.exp1.causal.model = df.exp1.counterfactual.long %>%
group_by(clip) %>%
summarize(empirical = mean(rating)) %>%
ungroup() %>%
left_join(df.exp1.counterfactual.model %>%
select(-c(noise, sse)),
by = "clip") %>%
rename(simulation = prediction) %>%
mutate(empirical = ifelse(outcome_actual == 1, 100 - empirical, empirical),
simulation = ifelse(outcome_actual == 1, 100 - simulation, simulation))`summarise()` ungrouping output (override with `.groups` argument)2.4 Counterfactual condition
2.4.1 Plots
set.seed(1)
df.plot = df.exp1.counterfactual.long %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.)/2)) %>%
left_join(df.exp1.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c(rating, prediction),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss", "actual close", "actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")),
model = factor(model, levels = c("rating", "prediction"))) %>%
mutate(clipindex = ifelse(clip == 11, 2, clipindex),
clipindex = ifelse(clip == 12, 1, clipindex)) #swap clips 11 and 12
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(1:18, each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
p = actual_counterfactual_plot(x = df.plot, ylabel = "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp1_counterfactual_bars.pdf",
# width = 8,
# height = 6)
2.4.2 Stats
2.4.2.1 Model evaluation
df.exp1.counterfactual.means %>%
left_join(df.exp1.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(noise = unique(noise),
simulation_r = cor(rating_mean, prediction),
simulation_rmse = rmse(rating_mean, prediction)) %>%
print_table()| noise | simulation_r | simulation_rmse |
|---|---|---|
| 1 | 0.97 | 10.95 |
2.5 Causal condition
2.5.1 Plots
set.seed(1)
model.name = c("empirical")
# model.name = c("simulation")
df.plot = df.exp1.causal.long %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp1.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c(rating, simulation, empirical),
names_to = "model",
values_to = "value") %>%
filter(model %in% c(model.name, "rating")) %>%
mutate(model = factor(model, levels = c("rating", model.name))) %>%
mutate(
colorindex = ifelse(outcome_actual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index.actual = factor(index_actual,
levels = c("actual miss", "actual close", "actual hit")),
index.counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit"))) %>%
mutate(clipindex = ifelse(clip == 11, 2, clipindex),
clipindex = ifelse(clip == 12, 1, clipindex)) # swap clips 11 and 12
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(1:18, each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
p = actual_counterfactual_plot(df.plot, "causal judgment")
print(p)
# ggsave("../../figures/plots/exp1_causal_bars.pdf",
# width = 8,
# height = 6)
2.5.2 Stats
2.5.2.1 Model evaluation
df.exp1.causal.long %>%
mutate(rating = abs(rating)) %>%
group_by(clip,
outcome_actual,
outcome_counterfactual,
index_actual,
index_counterfactual) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp1.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(empirical_r = cor(rating, empirical),
empirical_rmse = rmse(rating, empirical),
simulation_r = cor(rating, simulation),
simulation_rmse = rmse(rating, simulation)) %>%
print_table()`summarise()` regrouping output by 'clip', 'outcome_actual', 'outcome_counterfactual', 'index_actual' (override with `.groups` argument)| empirical_r | empirical_rmse | simulation_r | simulation_rmse |
|---|---|---|---|
| 0.96 | 8.57 | 0.93 | 15.15 |
2.5.2.2 Bayesian regression
2.5.2.3 Specific hypotheses
df.exp1.causal.posterior = df.exp1.causal.long %>%
group_by(clip, index_actual, index_counterfactual, outcome_actual) %>%
summarize(value = mean(rating)) %>%
ungroup() %>%
add_fitted_draws(newdata = .,
model = fit.brm_exp1_causal,
re_formula = NA) %>%
ungroup() %>%
select(clip, .value, .draw) %>%
pivot_wider(names_from = clip,
values_from = .value)`summarise()` regrouping output by 'clip', 'index_actual', 'index_counterfactual' (override with `.groups` argument)func_posterior_difference = function(data, clip1, clip2){
data %>%
mutate(difference = .data[[clip1]] - .data[[clip2]]) %>%
pull(difference) %>%
mean_hdci() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) %>%
print()
}
func_posterior_difference(data = df.exp1.causal.posterior,
clip1 = "8",
clip2 = "13") difference
1 (0.26 [-6.51, 7.31])2.5.3 Tables
fit.brm_exp1_causal %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(term = str_remove(term, "b_"),
term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
# print_table(format = "latex")
print_table()| term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|
| intercept | -18.33 | 3.36 | -23.81 | -12.78 |
| index_actualactualclose | 4.31 | 3.40 | -1.40 | 9.83 |
| index_actualactualhit | 4.06 | 4.84 | -4.01 | 11.88 |
| index_counterfactualcounterfactualclose | -17.52 | 3.37 | -22.96 | -12.01 |
| index_counterfactualcounterfactualhit | -61.37 | 4.38 | -68.54 | -54.30 |
| outcome_actual | 92.36 | 5.11 | 84.04 | 100.96 |
3 Experiment 2: Brick & Teleport
3.1 Clips
3.2 Read in data
df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv",
header = T,
stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv",
header = T,
stringsAsFactors = F)
# Information about each clip
df.exp2.clipinfo = tibble(clip = 1:20,
outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1,
0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0,
1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
index_actual = c(rep(c("actual miss",
"actual close",
"actual hit"),
each = 6),
rep("practice", 2)),
index_counterfactual = c(rep(rep(c("counterfactual miss",
"counterfactual close",
"counterfactual hit"),
each = 2),
3),
rep("practice", 2))) %>%
mutate(index_actual = factor(index_actual, levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
# Structure data
df.exp2.causal_first.long = df.exp2.causal_first$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:40, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(question = rep(rep(c("causal", "counterfactual"), each = 20), max(participant)),
condition = "causal_first") %>%
left_join(df.exp2.clipinfo,
by = "clip") %>%
arrange(participant, question, clip) %>%
select(participant, question, clip, rating, everything())
df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>%
enframe() %>%
mutate(participant = 1:n()) %>%
separate_rows(value, sep = ",") %>%
mutate(name = rep(c("clip", "rating"), n()/2),
order = rep(rep(1:40, each = 2), max(participant)),
value = as.numeric(value)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(question = rep(rep(c("counterfactual", "causal"), each = 20), max(participant)),
condition = "counterfactual_first") %>%
left_join(df.exp2.clipinfo,
by = "clip") %>%
arrange(participant, question, clip) %>%
select(participant, question, clip, rating, everything())
# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
bind_rows(df.exp2.counterfactual_first.long %>%
mutate(participant = participant +
max(df.exp2.causal_first.long$participant)))3.2.1 Demographics
# counterfactual condition
df.exp2.counterfactual_first %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
sep = ",") %>%
select(gender:time) %>%
bind_rows(df.exp2.causal_first %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
sep = ",") %>%
select(gender:time)) %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 82 | 35 | 12.2 | 45 | 21.25 | 5.11 |
3.3 Read in model predictions
3.3.1 Counterfactual condition
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "teleport")]
df.exp2.model = files %>%
map_dfr(~ read.csv(str_c(path, .))) %>%
set_names(c("clip", "noise", "prediction")) %>%
left_join(df.exp2.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))
# calculate mean counterfactual judgments
df.exp2.counterfactual.means = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
# mutate(rating = abs(rating)) %>%
group_by(clip) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp2.clipinfo, by = "clip")`summarise()` ungrouping output (override with `.groups` argument)# find noisy simulation model that best predicts the mean counterfactual judgments by
# minimizing the sum of squared errors
df.exp2.counterfactual.model = df.exp2.model %>%
group_by(noise) %>%
nest() %>%
mutate(sse = map(data,
~ sum((.x$prediction*100 - df.exp2.counterfactual.means$rating) ^ 2)),
sse = unlist(sse)) %>%
ungroup() %>%
filter(sse == min(sse)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)3.3.2 Causal condition
# Model predictions based on counterfactual judgments
df.exp2.causal.model = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
group_by(clip, condition) %>%
summarize(rating = mean(rating)) %>%
pivot_wider(names_from = condition,
values_from = rating) %>%
left_join(df.exp2.combined.long %>%
filter(question == "counterfactual", clip <= 18) %>%
group_by(clip) %>%
summarize(combined = mean(rating)),
by = "clip") %>%
left_join(df.exp2.counterfactual.model,
by = "clip") %>%
select(-c(sse, noise)) %>%
rename(simulation = prediction) %>%
mutate_at(.vars = vars(simulation, causal_first, counterfactual_first, combined),
.funs = list(~ ifelse(outcome_actual == 1, 100 - ., .)))`summarise()` regrouping output by 'clip' (override with `.groups` argument)`summarise()` ungrouping output (override with `.groups` argument)3.4 Counterfactual condition
3.4.1 Plots
set.seed(1)
df.plot = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
mutate(clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp2.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c(rating, prediction),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")),
model = factor(model,
levels = c("rating", "prediction")))
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
"7", "8 (C)", "9", "10", "11", "12 (C)",
"13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"), each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
p = actual_counterfactual_plot(df.plot, "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp2_counterfactual_bars.pdf",
# width = 8,
# height = 6)
3.4.2 Stats
3.4.2.1 Model evaluation
df.exp2.counterfactual.means %>%
left_join(df.exp2.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(noise = unique(noise),
simulation_r = cor(rating, prediction),
simulation_rmse = rmse(rating, prediction)) %>%
print_table()| noise | simulation_r | simulation_rmse |
|---|---|---|
| 0.9 | 0.96 | 13.46 |
3.4.2.2 Bayesian regressions
3.4.2.2.1 Check for potential order effects
# random intercepts & slopes model
fit.brm_exp2_counterfactual = brm(formula = rating ~ 1 + condition *
(index_actual + index_counterfactual) +
(1 + index_actual +
index_counterfactual | participant),
data = df.exp2.combined.long %>%
na.omit() %>%
filter(question == "counterfactual"),
cores = 2,
seed = 1,
file = "cache/brm_exp2_counterfactual")3.4.3 Tables
fit.brm_exp2_counterfactual %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(term = str_remove(term, "b_"),
term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
# print_table(format = "latex")
print_table()| term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|
| intercept | 13.38 | 2.58 | 9.12 | 17.55 |
| conditioncounterfactual_first | -1.11 | 3.65 | -6.99 | 4.82 |
| index_actualactualclose | -1.48 | 2.78 | -6.02 | 3.11 |
| index_actualactualhit | -4.42 | 2.40 | -8.32 | -0.46 |
| index_counterfactualcounterfactualclose | 37.94 | 3.33 | 32.39 | 43.35 |
| index_counterfactualcounterfactualhit | 73.80 | 4.17 | 67.00 | 80.57 |
| conditioncounterfactual_first:index_actualactualclose | -0.55 | 3.84 | -7.11 | 5.73 |
| conditioncounterfactual_first:index_actualactualhit | -0.85 | 3.41 | -6.52 | 4.85 |
| conditioncounterfactual_first:index_counterfactualcounterfactualclose | 1.23 | 4.79 | -6.62 | 9.18 |
| conditioncounterfactual_first:index_counterfactualcounterfactualhit | -0.02 | 5.92 | -9.86 | 9.55 |
3.5 Causal condition
3.5.1 Plots
set.seed(1)
func_exp2_causal_plot = function(condition.name, model.name){
df.plot = df.exp2.combined.long %>%
filter(question == "causal", clip <= 18) %>%
filter(condition %in% condition.name) %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp2.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
pivot_longer(cols = c("rating", all_of(model.name)),
names_to = "model",
values_to = "value") %>%
mutate(model = factor(model, levels = c("rating", model.name))) %>%
mutate(colorindex = ifelse(outcome_actual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
"7", "8 (C)", "9", "10", "11", "12 (C)",
"13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"), each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA)
actual_counterfactual_plot(df.plot, "causal judgment")
}
plot.list = list(func_exp2_causal_plot(condition.name = "counterfactual_first",
model.name = "combined"),
func_exp2_causal_plot(condition.name = "causal_first",
model.name = "combined"))
# model.name = "combined"
# condition.name = "counterfactual_first"
# condition.name = "causal_first"
#
# func_exp2_causal_plot(condition.name = condition.name,
# model.name = model.name),
#
# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
# width = 8,
# height = 6)3.5.2 Stats
3.5.2.1 Model evaluation
df.data = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18) %>%
mutate(rating = abs(rating)) %>%
group_by(condition,
clip,
outcome_actual,
outcome_counterfactual,
index_actual,
index_counterfactual) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp2.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual"))`summarise()` regrouping output by 'condition', 'clip', 'outcome_actual', 'outcome_counterfactual', 'index_actual' (override with `.groups` argument)df.data %>%
group_by(condition) %>%
summarize(empirical_r = cor(rating, combined),
empirical_rmse = rmse(rating, combined)) %>%
print_table()`summarise()` ungrouping output (override with `.groups` argument)| condition | empirical_r | empirical_rmse |
|---|---|---|
| causal_first | 0.92 | 13.78 |
| counterfactual_first | 0.99 | 5.27 |
df.data %>%
summarize(empirical_r = cor(rating, combined),
empirical_rmse = rmse(rating, combined),
simulation_r = cor(rating, simulation),
simulation_rmse = rmse(rating, simulation)) %>%
print_table()| empirical_r | empirical_rmse | simulation_r | simulation_rmse |
|---|---|---|---|
| 0.95 | 10.43 | 0.92 | 19.09 |
3.5.2.2 Bayesian regression
# random intercepts + slopes model
fit.brm_exp2_causal = brm(formula = rating ~ 1 + condition *
(index_actual + index_counterfactual + outcome_actual) +
(1 + index_actual + index_counterfactual +
outcome_actual | participant),
data = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18),
cores = 2,
seed = 1,
file = "cache/brm_exp2_causal")3.5.2.3 Specific hypotheses
df.exp2.causal.posterior = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18) %>%
group_by(clip, index_actual, index_counterfactual, outcome_actual, condition) %>%
summarize(value = mean(rating)) %>%
ungroup() %>%
add_fitted_draws(newdata = .,
model = fit.brm_exp2_causal,
re_formula = NA) %>%
ungroup() %>%
select(clip, condition, .value, .draw) %>%
pivot_wider(names_from = clip,
values_from = .value)`summarise()` regrouping output by 'clip', 'index_actual', 'index_counterfactual', 'outcome_actual' (override with `.groups` argument)# difference between clips 13 and 14 versus 8
df.exp2.causal.posterior %>%
mutate(difference = (`13` + `14`) / 2 - `8`) %>%
pull(difference) %>%
mean_hdi() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) difference
1 (-0.55 [-12.09, 10.97])# difference between clips 5 and 6 versus 11
df.exp2.causal.posterior %>%
mutate(difference = (`5` + `6`) / 2 - `11`) %>%
pull(difference) %>%
mean_hdi() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) difference
1 (-5.48 [-13.8, 2.73])3.5.3 Tables
3.5.3.1 Bayesian regression
fit.brm_exp2_causal %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(term = str_remove(term, "b_"),
term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
# print_table(format = "latex")
print_table()| term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|
| intercept | -24.07 | 3.50 | -29.80 | -18.27 |
| conditioncounterfactual_first | 11.47 | 4.91 | 3.29 | 19.40 |
| index_actualactualclose | 8.17 | 3.44 | 2.39 | 13.62 |
| index_actualactualhit | 13.42 | 4.90 | 5.42 | 21.30 |
| index_counterfactualcounterfactualclose | -30.20 | 3.72 | -36.41 | -24.36 |
| index_counterfactualcounterfactualhit | -50.29 | 4.60 | -57.94 | -42.84 |
| outcome_actual | 98.53 | 5.40 | 89.53 | 107.57 |
| conditioncounterfactual_first:index_actualactualclose | -5.39 | 4.94 | -13.38 | 2.80 |
| conditioncounterfactual_first:index_actualactualhit | -16.97 | 7.02 | -28.61 | -5.75 |
| conditioncounterfactual_first:index_counterfactualcounterfactualclose | -8.13 | 5.17 | -16.51 | 0.35 |
| conditioncounterfactual_first:index_counterfactualcounterfactualhit | -25.51 | 6.44 | -36.53 | -15.18 |
| conditioncounterfactual_first:outcome_actual | 10.72 | 7.60 | -1.63 | 23.34 |
4 Experiment 3: 3 balls
4.1 Clips
4.2 Read in data
# clipinfo
df.exp3.clipinfo = tibble(clip = 1:32,
outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1,
0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1,
0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
clipindex = rep(1:2, 16))
# COUNTERFACTUAL JUDGMENTS
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv",
header = T,
stringsAsFactors = F)
# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("gender", "age", "difficulty", "smoothness", "time", "condition"),
n()/6),
participant = rep(1:(n()/6), each = 6)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(condition = factor(condition, levels = 20:21, labels = c("A", "B")),
gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>%
select(participant, condition, gender, age, time, smoothness, difficulty)
# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("clip", "rating"), n()/2),
participant = rep(1:(n()/64), each = 64),
order = rep(rep(1:32, each = 2), n()/64)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
left_join(df.exp3.counterfactual.demographics %>%
select(participant, ball = condition),
by = "participant") %>%
select(participant, ball, clip, everything()) %>%
arrange(participant, clip)
# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv",
header = T,
stringsAsFactors = F)
# demographics
df.exp3.causal.demographics = df.exp3.causal$extra %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("gender", "age", "difficulty", "smoothness", "time",
"counterbalance", "replayType", "experiment"),
n()/8),
participant = rep(1:(n()/8), each = 8)) %>%
pivot_wider(names_from = name,
values_from = value) %>%
mutate(counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>%
select(participant, gender, age, counterbalance, time, smoothness, difficulty)
# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
enframe() %>%
separate_rows(value) %>%
mutate(value = as.numeric(value),
name = rep(c("clip", "x", "y", "z", "rating1", "rating2"), n()/6),
participant = rep(1:(n()/(32*6)), each = (32*6)),
order = rep(rep(1:32, each = 6), n()/(32*6))) %>%
filter(!name %in% c("x", "y", "z")) %>%
pivot_wider(names_from = name,
values_from = value) %>%
left_join(df.exp3.causal.demographics %>%
select(participant, counterbalance),
by = "participant") %>%
mutate(A = ifelse(counterbalance == 1, rating1, rating2),
B = ifelse(counterbalance == 1, rating2, rating1)) %>%
select(-rating1, -rating2) %>%
pivot_longer(cols = c(A, B),
names_to = "ball",
values_to = "rating") %>%
select(participant, clip, ball, rating, order) %>%
arrange(participant, clip, ball)4.2.1 Demographics
# counterfactual condition
df.exp3.counterfactual.demographics %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 80 | 33 | 10.1 | 34 | 18.08 | 4.63 |
# causal condition
df.exp3.causal.demographics %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) %>%
print_table()| n | age_mean | age_sd | n_female | time_mean | time_sd |
|---|---|---|---|---|---|
| 41 | 34 | 10.5 | 21 | 21.19 | 4.96 |
4.3 Read in model predictions
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "3ball")]
df.exp3.model = files %>%
map_dfr(~ read_csv(str_c(path, .))) %>%
rename(clip = trial)4.3.1 Counterfactual condition
# calculate mean counterfactual judgments
df.exp3.counterfactual.means = df.exp3.counterfactual.long %>%
group_by(clip, ball) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup()`summarise()` regrouping output by 'clip' (override with `.groups` argument)# find noisy simulation model that best predicts the mean counterfactual judgments by
# minimizing the sum of squared errors
df.exp3.model.counterfactual = df.exp3.model %>%
group_by(noise) %>%
select(clip, contains("whether"), noise) %>%
pivot_longer(cols = c(A_whether, B_whether),
names_to = "ball",
values_to = "prediction") %>%
mutate(ball = str_remove(ball, "_whether")) %>%
arrange(noise, clip, ball) %>%
left_join(df.exp3.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_both == 1, 1 - prediction, prediction)) %>%
group_by(noise) %>%
nest() %>%
mutate(sse = map(data,
~ sum((.x$prediction*100 -
df.exp3.counterfactual.means$rating_mean) ^ 2)),
sse = unlist(sse)) %>%
ungroup() %>%
filter(sse == min(sse)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)4.3.2 Causal condition
# calculate mean causal judgments
df.exp3.causal.means = df.exp3.causal.long %>%
group_by(clip, ball) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup()`summarise()` regrouping output by 'clip' (override with `.groups` argument)df.exp3.model.causal = df.exp3.model %>%
# take into account difference-making
mutate_at(.vars = vars(A_how:A_robust),
.funs = list(~ . * A_difference)) %>%
mutate_at(.vars = vars(B_how:B_robust),
.funs = list(~ . * B_difference)) %>%
# choose model based on best fit with counterfactual data
filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
pivot_longer(cols = A_difference:B_robust,
names_to = c("ball", "aspect"),
names_sep = "_",
values_to = "value") %>%
pivot_wider(names_from = aspect,
values_from = value) %>%
select(clip, ball:robust)4.3.3 Heuristic model
l.features = fromJSON("data/features.json")
df.features.initial = l.features[["appearance"]] %>%
cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
setNames(c(str_c("ball", LETTERS[c(1, 2, 5)], "_appearance"),
"ballA_velx",
"ballA_vely",
"ballB_velx",
"ballB_vely",
"ballE_velx",
"ballE_vely")) %>%
mutate(clip = 1:nrow(.)) %>%
pivot_longer(cols = starts_with("ball")) %>%
separate(name, into = c("ball", "property")) %>%
pivot_wider(names_from = property, values_from = value) %>%
mutate(ball = str_remove_all(ball, pattern = "ball"))
# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")
for (i in 1:32) {
df.tmp = data.frame()
for (j in 1:length(ballnames)) {
ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
if (ncollisions > 0) {
tmp = data.frame(
ball = ballnames[j],
object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
ncollision = 1:ncollisions
)
df.tmp = df.tmp %>%
rbind(tmp)
}
}
df.features.collisions = df.features.collisions %>%
rbind(df.tmp %>%
mutate(clip = i) %>%
select(clip, ball, everything()))
}
# find end of clip
tmp = df.features.collisions %>%
group_by(clip) %>%
filter(ball == "ballE",
str_detect(object, "goal|Left|Right")) %>%
group_by(clip) %>%
mutate(endclip = max(time)) %>%
select(clip, endclip) %>%
ungroup() %>%
distinct()
# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
left_join(tmp, by = "clip") %>%
mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
group_by(clip) %>%
filter(time <= endclip)
# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
filter(ball == "E") %>%
mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
select(clip, E.moving)
# contact with ball E
tmp.contactE = df.features.collisions %>%
filter(object == "ballE") %>%
mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
mutate(contactE = 1) %>%
select(clip, ball, contactE) %>%
group_by(clip, ball) %>%
summarize(contactE = any(contactE %>% as.logical()) * 1) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
by = c("clip", "ball")) %>%
mutate(contactE = ifelse(is.na(contactE), 0, contactE))
# number of collisions
tmp.ncollisions = df.features.collisions %>%
filter(str_detect(object, "ball"),
ball != "ballE") %>%
mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
group_by(clip, ball) %>%
count() %>%
rename(ncollision = n)
# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
filter(str_detect(object, "ball"),
ball != "ballE") %>%
count(clip, ball, object) %>%
mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)) %>%
group_by(clip) %>%
summarize(AE = any(AE %>% as.logical()) * 1,
BE = any(BE %>% as.logical()) * 1) %>%
mutate(A = ifelse(AE == 1 & BE == 0, 1, 0),
B = ifelse(AE == 0 & BE == 1, 1, 0)) %>%
select(clip, A, B) %>%
pivot_longer(cols = -clip,
names_to = "ball",
values_to = "exclusive")
# impact
func_angle = function(x, y) {
dot.prod = x %*% y
norm.x = norm(x, type = "2")
norm.y = norm(y, type = "2")
theta = acos(dot.prod / (norm.x * norm.y))
as.numeric(theta)
}
# impact on ball E
tmp.impactE = df.features.collisions %>%
rowwise() %>%
filter(str_detect(object, "ball"),
ball == "ballE") %>%
mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
post.speed = abs(post.velx) + abs(post.vely),
speed.diff = post.speed - pre.speed,
direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
direction.diff = ifelse(is.na(direction.diff),
abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
direction.diff)) %>%
ungroup() %>%
group_by(clip, object) %>%
summarize(speed.diff = sum(speed.diff),
direction.diff = sum(direction.diff)) %>%
rename(ball = object) %>%
mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
by = c("clip", "ball")) %>%
mutate_at(.vars = vars(speed.diff, direction.diff),
.funs = ~ ifelse(is.na(.), 0, .)) %>%
arrange(clip, ball) %>%
rename(E.speed.diff = speed.diff,
E.direction.diff = direction.diff)
# total impact
tmp.impactTotal = df.features.collisions %>%
rowwise() %>%
filter(str_detect(object, "ballA|ballB")) %>%
mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
post.speed = abs(post.velx) + abs(post.vely),
speed.diff = post.speed - pre.speed,
direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
direction.diff = ifelse(is.na(direction.diff),
abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
direction.diff)) %>%
ungroup() %>%
group_by(clip, object) %>%
summarize(speed.diff = sum(speed.diff),
direction.diff = sum(direction.diff)) %>%
rename(ball = object) %>%
mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
by = c("clip", "ball")) %>%
mutate_at(.vars = vars(speed.diff, direction.diff),
.funs = ~ ifelse(is.na(.), 0, .)) %>%
arrange(clip, ball) %>%
rename(total.speed.diff = speed.diff,
total.direction.diff = direction.diff)
# transfer of force
tmp.transfer = df.features.collisions %>%
filter(str_detect(ball, "ballA|ballB"),
str_detect(object, "ball")) %>%
mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)) %>%
mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
arrange(clip, time) %>%
group_by(clip) %>%
summarize(A = any(!is.na(AE)),
B = any(!is.na(BE)),
A = ifelse(any(!is.na(AB)) &
any(!is.na(BE)) &
max(AB, na.rm = T) < min(BE, na.rm = T),
T,
A),
B = ifelse(any(!is.na(AB)) &
any(!is.na(AE)) &
max(AB, na.rm = T) < min(AE, na.rm = T),
T,
B)) %>%
pivot_longer(cols = -clip,
names_to = "ball",
values_to = "transfer") %>%
arrange(clip, ball)
# collect features
df.features = df.features.initial %>%
filter(ball != "E") %>%
mutate(ball = factor(ball)) %>%
mutate(moving = ifelse(velx == 0 & vely == 0, 0, 1),
speed = abs(velx) + abs(vely)) %>%
left_join(tmp.contactE) %>%
left_join(tmp.impactE) %>%
left_join(tmp.impactTotal) %>%
left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
left_join(tmp.movingE) %>%
left_join(tmp.exclusive) %>%
select(-c(appearance, velx, vely))
# remove temporary variables
rm(list = ls()[str_detect(ls(), "tmp")])4.4 Counterfactual condition
4.4.1 Plots
4.4.1.1 Bar graph
set.seed(1)
df.plot = df.exp3.counterfactual.long %>%
left_join(df.exp3.model.counterfactual %>%
select(clip, ball, prediction) %>%
mutate(ball = as.factor(ball)),
by = c("clip", "ball")) %>%
left_join(df.exp3.clipinfo,
by = "clip") %>%
pivot_longer(cols = c(rating, prediction),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = 1,
colorindex = ifelse(model == "rating" & ball == "A" & outcome_b == 1,
2,
colorindex),
colorindex = ifelse(model == "rating" & ball == "B" & outcome_a == 1,
2,
colorindex),
colorindex = ifelse(model == "prediction", 3, colorindex),
colorindex = as.factor(colorindex),
model = factor(model, levels = c("rating", "prediction"))) %>%
arrange(participant, clip, ball)
df.text = df.plot %>%
expand(ball, clip, model) %>%
mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
y = 110,
x = 1.5,
colorindex = NA)
p = actual_counterfactual_threeballs_plot(df.plot, "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp3_counterfactual_bars.pdf",
# width = 12,
# height = 6)
4.4.2 Stats
4.4.2.1 Model evaluation
# best fitting model
df.exp3.counterfactual.means %>%
left_join(df.exp3.model.counterfactual,
by = c("clip", "ball")) %>%
summarize(simulation_r = cor(rating_mean, prediction),
simulation_rmse = rmse(rating_mean, prediction)) %>%
print_table()| simulation_r | simulation_rmse |
|---|---|
| 0.86 | 19.84 |
# deterministic model
df.exp3.counterfactual.means %>%
left_join(df.exp3.clipinfo %>%
select(clip, outcome_a, outcome_b) %>%
pivot_longer(cols = -clip,
names_to = "ball",
values_to = "outcome") %>%
mutate(ball = factor(ball,
levels = c("outcome_a", "outcome_b"),
labels = c("B", "A"))),
by = c("clip", "ball")) %>%
mutate(outcome = outcome * 100) %>%
summarize(simulation_r = cor(rating_mean, outcome),
simulation_rmse = rmse(rating_mean, outcome)) %>%
print_table()| simulation_r | simulation_rmse |
|---|---|
| 0.83 | 30.28 |
4.5 Causal condition
4.5.1 Stats
4.5.1.1 Bayesian mixed effects model
df.data = df.exp3.causal.long %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball"))
# using whether-causation from the model
fit.brm_exp3_causal_w = brm(formula = rating ~ 1 + whether + (1 + whether | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_w")
fit.brm_exp3_causal_wh = brm(formula = rating ~ 1 + whether + how +
(1 + whether + how | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_wh")
fit.brm_exp3_causal_whs = brm(formula = rating ~ 1 + whether + how + sufficient +
(1 + whether + how + sufficient | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_whs")
# fit.brm_exp3_causal_w %>% summary()
# fit.brm_exp3_causal_wh %>% summary()
# fit.brm_exp3_causal_whs %>% summary()4.5.1.2 Model comparison
fit.brm_exp3_causal_w = add_criterion(fit.brm_exp3_causal_w,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.fit.brm_exp3_causal_wh = add_criterion(fit.brm_exp3_causal_wh,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.fit.brm_exp3_causal_whs = add_criterion(fit.brm_exp3_causal_whs,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object. elpd_diff se_diff
fit.brm_exp3_causal_whs 0.0 0.0
fit.brm_exp3_causal_wh -125.2 15.6
fit.brm_exp3_causal_w -426.2 31.1 4.5.1.3 Heuristic model
- z-scored predictors and invidiual responses
# Regression based on features
df.regression.features = df.exp3.causal.long %>%
left_join(df.exp3.clipinfo, by = "clip") %>%
left_join(df.features %>%
mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
by = c("clip", "ball")) %>%
clean_names() %>%
mutate(e_moving = 1 - e_moving)
# restrict the weights such that all predictors are positive
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))
fit.brm_exp3_causal_heuristic = brm(formula = rating ~ moving +
speed +
contact_e +
e_speed_diff +
e_direction_diff +
total_speed_diff +
total_direction_diff +
transfer +
e_moving +
exclusive + (1 | participant),
prior = priors,
data = df.regression.features %>%
select(-c(clip, ball, order, clipindex,
contains("outcome"))),
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_heuristic")
fit.brm_exp3_causal_heuristic Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant)
Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 41)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 8.53 1.22 6.44 11.18 1.00 1163 1936
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 49.79 1.46 46.91 52.74 1.00 1129
moving 0.22 0.22 0.01 0.81 1.00 3399
speed 2.06 0.85 0.44 3.67 1.00 2319
contact_e 0.38 0.35 0.01 1.29 1.00 3201
e_speed_diff 0.12 0.12 0.00 0.46 1.00 4435
e_direction_diff 1.04 0.74 0.04 2.71 1.00 2202
total_speed_diff 2.21 0.95 0.41 4.17 1.00 2208
total_direction_diff 3.99 0.91 2.14 5.68 1.00 2718
transfer 15.59 0.79 14.06 17.14 1.00 3555
e_moving 0.18 0.17 0.01 0.62 1.00 3296
exclusive 4.39 0.73 2.97 5.81 1.00 3526
Tail_ESS
Intercept 1854
moving 1547
speed 1266
contact_e 2020
e_speed_diff 2268
e_direction_diff 1685
total_speed_diff 1393
total_direction_diff 1742
transfer 2992
e_moving 1738
exclusive 2538
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 33.21 0.46 32.33 34.10 1.00 5313 2845
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).4.5.1.4 Predictions for mean causal judgments
func_fit_data = function(fit){
fit %>%
fitted(newdata = df.exp3.model.causal %>%
left_join(df.features %>%
mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
by = c("clip", "ball")) %>%
clean_names() %>%
mutate(e_moving = 1 - e_moving),
re_formula = NA) %>%
as_tibble() %>%
pull(Estimate)
}
df.exp3.causal.regression = df.exp3.causal.means %>%
mutate(w = func_fit_data(fit.brm_exp3_causal_w),
wh = func_fit_data(fit.brm_exp3_causal_wh),
whs = func_fit_data(fit.brm_exp3_causal_whs),
heuristic = func_fit_data(fit.brm_exp3_causal_heuristic))4.5.1.5 Model evaluation
prediction = "whs"
df.exp3.causal.regression %>%
summarize(simulation_r = cor(rating_mean, .[[prediction]]),
simulation_rmse = rmse(rating_mean, .[[prediction]])) %>%
print_table()| simulation_r | simulation_rmse |
|---|---|
| 0.87 | 13.07 |
4.5.1.6 Individual participant regressions
df.fit = df.exp3.causal.long %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball"))
# priors
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))
# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline =
brm(formula = rating ~ 1,
data = df.fit %>%
filter(participant == 1),
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_baseline"))
fit.brm_exp3_causal_individual_w =
brm(formula = rating ~ 1 + whether,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_w"))
fit.brm_exp3_causal_individual_wh =
brm(formula = rating ~ 1 + whether + how,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_wh"))
fit.brm_exp3_causal_individual_whs =
brm(formula = rating ~ 1 + whether + how + sufficient,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_whs"))
# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
group_by(participant) %>%
nest() %>%
ungroup() %>%
# fit model for each participant
mutate(fit_baseline = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_baseline,
newdata = .x)),
fit_w = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_w,
newdata = .x)),
fit_wh = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_wh,
newdata = .x)),
fit_whs = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_whs,
newdata = .x))) %>%
mutate(fit_baseline = map(.x = fit_baseline,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_w = map(.x = fit_w,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_wh = map(.x = fit_wh,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_whs = map(.x = fit_whs,
~ add_criterion(.x, criterion = c("loo", "waic"))),
r_w = map2_dbl(.x = data,
.y = fit_w,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
r_wh = map2_dbl(.x = data,
.y = fit_wh,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
r_whs = map2_dbl(.x = data,
.y = fit_whs,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_w = map2_dbl(.x = data,
.y = fit_w,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_wh = map2_dbl(.x = data,
.y = fit_wh,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_whs = map2_dbl(.x = data,
.y = fit_whs,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
model_comparison = pmap(.l = list(baseline = fit_baseline,
w = fit_w,
wh = fit_wh,
whs = fit_whs),
.f = ~ loo_compare(..1, ..2, ..3, ..4)),
best_model = map_chr(.x = model_comparison,
.f = ~ rownames(.) %>% .[1]),
best_model = factor(best_model,
levels = c("..1", "..2", "..3", "..4"),
labels = c("baseline", "w", "wh", "whs")))
save(list = c("df.exp3.causal.individual_fit"),
file = "data/exp3_causal_individual_fit.RData")4.5.1.7 Load individual participant regressions
load("data/exp3_causal_individual_fit.RData")
# count how many participants are best fit by the different models
df.exp3.causal.individual_fit %>%
count(best_model) %>%
print_table()| best_model | n |
|---|---|
| wh | 2 |
| whs | 39 |
4.5.1.8 Clusters of individual participants’ responses
set.seed(1)
df.cluster_whs = fit.brm_exp3_causal_whs %>%
ranef() %>%
.$participant %>%
as_tibble() %>%
select(contains("Estimate")) %>%
set_names(tolower(str_replace_all(names(.), "Estimate.", ""))) %>%
mutate(participant = 1:n()) %>%
relocate(participant)
df.cluster_whs = df.cluster_whs %>%
mutate(cluster = df.cluster_whs %>%
select(-participant) %>%
kmeans(centers = 2) %>%
.$cluster)4.5.2 Plots
4.5.2.1 Bar graph
set.seed(1)
model_index = "whs"
# model predictions
model_prediction = list(fit.brm_exp3_causal_w,
fit.brm_exp3_causal_wh,
fit.brm_exp3_causal_whs) %>%
map_dfr(~ fitted(., df.exp3.causal.means %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball")),
re_formula = NA) %>%
as_tibble()) %>%
mutate(ball = rep(c("A", "B"), n()/2),
clip = rep(rep(1:32, each = 2), 3),
model = rep(c("w", "wh", "whs"),
each = 64))
df.plot = df.exp3.causal.long %>%
left_join(df.exp3.clipinfo %>%
select(clip, outcome_both),
by = c("clip")) %>%
left_join(model_prediction %>%
filter(model == model_index) %>%
select(model = Estimate, clip, ball),
by = c("clip", "ball")) %>%
pivot_longer(cols = c(rating, model),
names_to = "model",
values_to = "value") %>%
mutate(colorindex = 1,
colorindex = ifelse(model == "rating" & outcome_both == 0, 1, colorindex),
colorindex = ifelse(model == "rating" & outcome_both == 1, 2, colorindex),
colorindex = ifelse(model == "model", 3, colorindex),
colorindex = as.factor(colorindex),
model = factor(model, levels = c("rating", "model"))) %>%
arrange(participant, clip, ball)
df.text = df.plot %>%
expand(ball, clip, model) %>%
mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
y = 110,
x = 1.5,
colorindex = NA)
p = actual_counterfactual_threeballs_plot(df.plot,"causal responsibility")
print(p)
# ggsave("../../figures/plots/exp3_causal_bars.pdf",
# width = 12,
# height = 6)
4.5.2.2 Selection
# \u2713 = check mark
# \u2717 = cross mark
check = "\u2713"
cross = "\u2717"
x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
df.plot = df.exp3.causal.long %>%
filter(clip %in% c(3, 7, 15, 16, 23)) %>%
group_by(clip, ball) %>%
summarize(mean = mean(rating),
low = smean.cl.boot(rating)[2],
high = smean.cl.boot(rating)[3]) %>%
left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
by = c("clip", "ball")) %>%
ungroup() %>%
pivot_longer(cols = c(mean, w, wh, whs),
names_to = "index",
values_to = "value") %>%
mutate_at(.vars = vars(low, high), .funs = ~ ifelse(index == "mean", ., NA)) %>%
mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
clip = factor(clip, levels = c(7, 23, 3, 15, 16),
labels = c("causal chain", "double prevention", "joint causation",
"overdetermination", "preemption")))`summarise()` regrouping output by 'clip' (override with `.groups` argument)df.labels = df.plot %>%
distinct(clip, ball) %>%
arrange(clip, ball) %>%
mutate(labels = x_labels,
value = -10,
index = NA)
func_load_image = function(clip){
readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}
df.clips = df.plot %>%
distinct(clip) %>%
arrange(clip) %>%
mutate(number = c(7, 23, 3, 15, 16),
grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
index = NA,
value = NA,
ball = NA,
label = str_c("clip ", number))
ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)
ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)
p = ggplot(data = df.plot,
mapping = aes(x = ball,
y = value,
group = index,
fill = index)) +
geom_bar(stat = "identity", color = "black",
position = position_dodge(0.9),
width = 0.9) +
geom_errorbar(mapping = aes(ymin = low, ymax = high),
width = 0,
size = 1,
position = position_dodge(0.9)) +
annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -25, ymax = -7) +
annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -25, ymax = -7) +
geom_text(data = df.labels,
mapping = aes(label = labels,
y = -Inf),
vjust = 1.2,
family = "Arial Unicode MS",
size = 8) +
geom_custom(data = df.clips,
mapping = aes(data = grob, x = 1.5, y = Inf),
grob_fun = function(x) rasterGrob(x,
interpolate = T,
vjust = -0.25)) +
geom_text(data = df.clips,
mapping = aes(label = label,
y = Inf,
x = -Inf),
size = 9,
hjust = -0.2,
vjust = -4) +
facet_grid(~clip) +
scale_fill_grey(start = 1,
end = 0,
labels = c("mean rating ", expression(CSM[W] ~ " "),
expression(CSM[WH] ~ " "),
expression(CSM[WHS] ~ " "))) +
labs(y = "causal responsibility", fill = "") +
coord_cartesian(clip = "off") +
theme_bw() +
theme(legend.text.align = 0,
text = element_text(size = 20),
panel.grid = element_blank(),
legend.position = "bottom",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
legend.text = element_text(size = 30),
strip.text = element_text(size = 30),
legend.spacing = unit(0.5, "cm"),
legend.background = element_rect(fill = "transparent"),
legend.margin = margin(t = 5, unit = "cm"),
plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))
print(p)
# ggsave("../../figures/plots/exp3_selection_bars.pdf",
# width = 20,
# height = 10,
# device = cairo_pdf)
4.5.2.3 Scatter plot
func_scatterplot = function(model){
if(model == "heuristic"){
xlabel = "Heuristic"
}else{
xlabel = bquote(CSM[.(toupper(model))])
}
l.models = list(w = fit.brm_exp3_causal_w,
wh = fit.brm_exp3_causal_wh,
whs = fit.brm_exp3_causal_whs,
heuristic = fit.brm_exp3_causal_heuristic)
df.data = df.exp3.causal.means %>%
left_join(df.exp3.model.causal, by = c("clip", "ball")) %>%
left_join(df.regression.features, by = c("clip", "ball"))
df.plot = l.models[[model]] %>%
fitted(newdata = df.data,
re_formula = NA) %>%
as_tibble() %>%
clean_names() %>%
bind_cols(df.data)
p = ggplot(data = df.plot,
mapping = aes(x = estimate,
y = rating_mean)) +
geom_abline(intercept = 0,
slope = 1,
linetype = 2) +
geom_smooth(aes(y = estimate,
ymin = q2_5,
ymax = q97_5),
stat = "identity",
color = "black") +
geom_linerange(size = 0.75,
mapping = aes(ymin = rating_low,
ymax = rating_high),
color = "gray50") +
geom_point(size = 2) +
scale_color_manual(values = c("black", "#e41a1c"),
guide = F) +
coord_fixed(xlim = c(0, 100),
ylim = c(0, 100)) +
labs(x = xlabel,
y = "responsibility judgment") +
annotate(geom = "text",
label = str_c(
"r = ", cor(df.plot$estimate, df.plot$rating_mean) %>%
round(2) %>%
as.character() %>%
str_sub(start = 2),
"\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>%
round(2)),
x = 5,
y = 95,
size = 6,
hjust = 0) +
scale_x_continuous(breaks = seq(0, 100, 25)) +
scale_y_continuous(breaks = seq(0, 100, 25)) +
theme(text = element_text(size = 20))
return(p)
}
# for creating and saving an individual scatter plot
# model = "w"
# model = "wh"
# model = "whs"
# model = "heuristic"
# func_scatterplot(model)
# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
# width = 5,
# height = 5)
list(func_scatterplot(model = "w"),
func_scatterplot(model = "wh"),
func_scatterplot(model = "whs"),
func_scatterplot(model = "heuristic")) %>%
wrap_plots(ncol = 2)
4.5.2.4 Individual participant variance for a selection of clips (clustered)
# \u2713 = check mark
# \u2717 = cross mark
check = "\u2713"
cross = "\u2717"
x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
df.plot = df.exp3.causal.long %>%
filter(clip %in% c(7, 23, 3, 15, 16)) %>%
mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
labels = c("causal chain", "double prevention", "joint causation",
"overdetermination", "preemption")))
df.plot = df.plot %>%
left_join(df.cluster_whs %>%
group_by(cluster) %>%
mutate(group = factor(cluster, levels = 1:2, labels = str_c("n = ", group_size(.)))) %>%
ungroup() %>%
select(participant, cluster, group),
by = "participant")
df.labels = df.plot %>%
distinct(clip, ball) %>%
arrange(clip, ball) %>%
mutate(labels = x_labels,
rating = -10,
group = NA,
participant = NA)
func_load_image = function(clip){
readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}
df.clips = df.plot %>%
distinct(clip) %>%
arrange(clip) %>%
mutate(number = c(7, 23, 3, 15, 16),
grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
group = NA,
participant = NA,
ball = NA,
label = str_c("clip ", number))
ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)
ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)
p = ggplot(data = df.plot,
mapping = aes(x = ball,
y = rating,
group = participant,
shape = group)) +
geom_line(mapping = aes(linetype = "individual",
color = group),
size = 1,
alpha = 0.3) +
geom_point(mapping = aes(color = group),
size = 1,
alpha = 0.3) +
stat_summary(mapping = aes(group = group,
color = group,
linetype = "mean"),
fun = "mean",
geom = "line",
size = 1.5) +
stat_summary(data = df.plot %>% filter(cluster == 1),
mapping = aes(group = group,
color = group),
fun.data = "mean_cl_boot",
geom = "pointrange",
size = 1,
shape = 19) +
stat_summary(data = df.plot %>% filter(cluster == 2),
mapping = aes(group = group,
color = group),
fun.data = "mean_cl_boot",
geom = "pointrange",
size = 1,
shape = 19) +
annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -30, ymax = -10) +
annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -30, ymax = -10) +
geom_text(data = df.labels,
mapping = aes(label = labels,
y = -Inf),
vjust = 1.2,
family = "Arial Unicode MS",
size = 8) +
geom_custom(data = df.clips,
mapping = aes(data = grob, x = 1.5, y = Inf),
grob_fun = function(x) rasterGrob(x,
interpolate = T,
vjust = -0.25)) +
geom_text(data = df.clips,
mapping = aes(label = label,
y = Inf,
x = -Inf),
size = 9,
hjust = -0.2,
vjust = -4) +
facet_wrap(~ clip,
ncol = 8) +
coord_cartesian(xlim = c(0.9, 2.1),
ylim = c(-5, 105),
expand = F,
clip = "off") +
scale_y_continuous(breaks = seq(0, 100, 25),
labels = seq(0, 100, 25)) +
scale_color_brewer(palette = "Set1",
guide = "none") +
labs(y = "causal responsibility rating",
linetype = "legend",
shape = "") +
theme_bw() +
guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
keywidth = unit(1.2, "cm")),
shape = guide_legend(override.aes = list(color = c("#377eb8", "#e41a1c"),
shape = c(19, 19),
alpha = c(1, 1),
size = c(5, 5)))) +
theme(panel.grid = element_blank(),
text = element_text(size = 20),
legend.position = c(0.505, 0.23),
axis.title.x = element_blank(),
legend.text = element_text(size = 20),
legend.title = element_text(size = 24),
legend.background = element_rect(fill = "transparent"),
strip.text = element_text(size = 30),
axis.text.x = element_blank(),
legend.key = element_rect(fill = "transparent"),
legend.box = "horizontal",
legend.spacing.x = unit(0.1, "cm"),
panel.spacing.x = unit(1, "cm"),
plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))
print(p)
# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
# plot = p,
# width = 20,
# height = 8.5,
# device = cairo_pdf)
4.5.3 Tables
4.5.3.1 Relationship between predictors
df.exp3.model %>%
filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
pivot_longer(cols = A_difference:B_robust,
names_to = c("ball", "aspect"),
names_sep = "_",
values_to = "value") %>%
pivot_wider(names_from = aspect,
values_from = value) %>%
select(difference, whether, how, sufficient, robust) %>%
correlate() %>%
shave() %>%
fashion() %>%
# print_table(format = "latex")
print_table()| rowname | difference | whether | how | sufficient | robust |
|---|---|---|---|---|---|
| difference | |||||
| whether | .50 | ||||
| how | .79 | .27 | |||
| sufficient | .21 | .10 | .36 | ||
| robust | .43 | .93 | .24 | .20 |
4.5.3.2 Population level predictors in the mixed effects models
fit.brm_exp3_causal_w %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(model = "CSM_w") %>%
bind_rows(fit.brm_exp3_causal_wh %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(model = "CSM_wh")) %>%
bind_rows(fit.brm_exp3_causal_whs %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(model = "CSM_whs")) %>%
mutate(term = tolower(term),
term = str_remove_all(term, "b_")) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
select(model, everything()) %>%
# print_table(format = "latex")
print_table()| model | term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|---|
| CSM_w | intercept | 31.68 | 1.61 | 29.02 | 34.37 |
| CSM_w | whether | 39.30 | 2.16 | 35.76 | 42.82 |
| CSM_wh | intercept | 10.29 | 1.56 | 7.80 | 12.83 |
| CSM_wh | whether | 28.40 | 2.68 | 23.96 | 32.70 |
| CSM_wh | how | 36.07 | 2.60 | 31.80 | 40.27 |
| CSM_whs | intercept | 10.43 | 1.59 | 7.82 | 13.01 |
| CSM_whs | whether | 27.93 | 2.60 | 23.67 | 32.18 |
| CSM_whs | how | 26.47 | 2.91 | 21.76 | 31.35 |
| CSM_whs | sufficient | 30.90 | 3.06 | 25.77 | 35.96 |
4.5.3.3 Individual participants regression results
df.exp3.causal.individual_fit %>%
select_if(~ (is.numeric(.) | is.factor(.))) %>%
mutate_at(.vars = vars(contains("_w")), .funs = ~ round(., 2)) %>%
select(participant, everything(), best_model) %>%
# print_table(format = "latex")
print_table()| participant | r_w | r_wh | r_whs | rmse_w | rmse_wh | rmse_whs | best_model |
|---|---|---|---|---|---|---|---|
| 1 | 0.29 | 0.38 | 0.48 | 30.00 | 29.03 | 27.74 | whs |
| 2 | 0.20 | 0.77 | 0.77 | 39.32 | 27.94 | 27.55 | whs |
| 3 | 0.37 | 0.78 | 0.79 | 38.90 | 28.00 | 26.95 | whs |
| 4 | 0.41 | 0.55 | 0.63 | 43.62 | 40.60 | 38.31 | whs |
| 5 | 0.45 | 0.53 | 0.54 | 39.23 | 37.32 | 36.81 | whs |
| 6 | 0.18 | 0.54 | 0.54 | 40.84 | 36.32 | 36.07 | whs |
| 7 | 0.49 | 0.61 | 0.64 | 32.96 | 29.93 | 29.15 | whs |
| 8 | 0.39 | 0.63 | 0.68 | 38.37 | 33.26 | 31.64 | whs |
| 9 | 0.40 | 0.72 | 0.76 | 35.86 | 28.34 | 26.39 | whs |
| 10 | 0.28 | 0.52 | 0.53 | 29.13 | 26.28 | 25.85 | whs |
| 11 | 0.51 | 0.56 | 0.60 | 33.66 | 32.19 | 31.17 | whs |
| 12 | 0.56 | 0.57 | 0.70 | 32.78 | 32.08 | 28.84 | whs |
| 13 | 0.56 | 0.71 | 0.74 | 29.83 | 25.40 | 24.21 | whs |
| 14 | 0.43 | 0.47 | 0.50 | 33.72 | 32.82 | 32.23 | whs |
| 15 | 0.19 | 0.35 | 0.37 | 39.77 | 38.29 | 37.94 | whs |
| 16 | 0.62 | 0.67 | 0.76 | 40.49 | 37.98 | 34.39 | whs |
| 17 | 0.42 | 0.66 | 0.71 | 28.57 | 23.82 | 22.42 | whs |
| 18 | 0.31 | 0.35 | 0.39 | 37.78 | 37.09 | 36.57 | whs |
| 19 | 0.37 | 0.60 | 0.62 | 33.41 | 29.23 | 28.62 | whs |
| 20 | 0.57 | 0.61 | 0.68 | 36.47 | 35.00 | 32.71 | whs |
| 21 | 0.43 | 0.50 | 0.57 | 30.95 | 29.61 | 28.32 | whs |
| 22 | 0.51 | 0.67 | 0.71 | 33.72 | 29.36 | 28.00 | whs |
| 23 | 0.33 | 0.45 | 0.52 | 38.76 | 37.00 | 35.57 | whs |
| 24 | 0.46 | 0.63 | 0.69 | 34.01 | 29.96 | 28.22 | whs |
| 25 | 0.61 | 0.66 | 0.70 | 30.60 | 28.91 | 27.55 | whs |
| 26 | 0.57 | 0.71 | 0.76 | 40.43 | 35.03 | 32.80 | whs |
| 27 | 0.46 | 0.52 | 0.66 | 43.12 | 41.49 | 38.12 | whs |
| 28 | 0.30 | 0.45 | 0.46 | 29.33 | 27.62 | 27.38 | whs |
| 29 | 0.48 | 0.60 | 0.65 | 38.93 | 36.03 | 34.48 | whs |
| 30 | 0.18 | 0.70 | 0.70 | 33.06 | 24.95 | 24.71 | whs |
| 31 | 0.48 | 0.59 | 0.69 | 37.54 | 34.72 | 31.76 | whs |
| 32 | 0.37 | 0.61 | 0.70 | 40.62 | 35.94 | 32.92 | whs |
| 33 | 0.50 | 0.56 | 0.64 | 38.64 | 36.79 | 34.85 | whs |
| 34 | 0.24 | 0.51 | 0.52 | 44.50 | 40.86 | 40.19 | whs |
| 35 | 0.15 | 0.32 | 0.30 | 34.56 | 33.36 | 33.37 | wh |
| 36 | 0.19 | 0.84 | 0.82 | 45.38 | 29.58 | 29.70 | wh |
| 37 | 0.55 | 0.63 | 0.67 | 36.04 | 33.36 | 32.11 | whs |
| 38 | 0.52 | 0.61 | 0.67 | 36.03 | 33.39 | 31.53 | whs |
| 39 | 0.47 | 0.77 | 0.77 | 29.61 | 21.75 | 21.58 | whs |
| 40 | 0.46 | 0.78 | 0.79 | 32.41 | 23.60 | 22.69 | whs |
| 41 | 0.32 | 0.50 | 0.58 | 32.90 | 30.38 | 28.69 | whs |
4.5.3.4 CSMwhs predictions for selection of cases
df.exp3.model.causal %>%
left_join(df.exp3.causal.regression,
by = c("clip", "ball")) %>%
filter(clip %in% c(7, 23, 16, 3, 15)) %>%
mutate_at(.vars = vars(difference, robust, sufficient, whether, heuristic),
.funs = ~ round(., 2)) %>%
select(clip, ball, difference, whether, how, sufficient, robust) %>%
mutate(clip = factor(clip, levels = c(7, 23, 16, 3, 15))) %>%
mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
.funs = ~ ifelse(. < 0.5,
str_c("xmark (", ., ")"),
str_c("cmark (", ., ")"))) %>%
arrange(clip) %>%
# print_table(format = "latex")
print_table()| clip | ball | difference | whether | how | sufficient | robust |
|---|---|---|---|---|---|---|
| 7 | A | cmark (1) | xmark (0.34) | cmark (1) | xmark (0) | xmark (0.25) |
| 7 | B | cmark (1) | cmark (1) | cmark (1) | cmark (0.67) | cmark (0.6) |
| 23 | A | xmark (0.05) | xmark (0) | xmark (0) | xmark (0) | xmark (0) |
| 23 | B | cmark (0.91) | cmark (0.79) | xmark (0) | xmark (0) | cmark (0.72) |
| 16 | A | cmark (1) | xmark (0.23) | cmark (1) | cmark (1) | xmark (0.35) |
| 16 | B | xmark (0) | xmark (0) | xmark (0) | xmark (0) | xmark (0) |
| 3 | A | cmark (1) | cmark (0.88) | cmark (1) | xmark (0.12) | cmark (0.76) |
| 3 | B | cmark (1) | cmark (0.89) | cmark (1) | xmark (0.11) | cmark (0.75) |
| 15 | A | cmark (1) | xmark (0.01) | cmark (1) | cmark (0.99) | xmark (0.1) |
| 15 | B | cmark (1) | xmark (0.01) | cmark (1) | cmark (1) | xmark (0.1) |
4.5.3.5 CSMwhs predictions for all cases
df.exp3.causal.regression %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball")) %>%
left_join(df.exp3.clipinfo %>%
select(-clipindex),
by = c("clip")) %>%
mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
.funs = ~ . * 100) %>%
mutate_at(.vars = vars(-ball), .funs = ~ round(., 0)) %>%
select(clip, ball, contains("outcome"), difference, whether, how, sufficient, robust,
w, wh, whs, heuristic, rating = rating_mean) %>%
# print_table(format = "latex", digits = 0)
print_table(digits = 0)| clip | ball | outcome_both | outcome_a | outcome_b | outcome_none | difference | whether | how | sufficient | robust | w | wh | whs | heuristic | rating |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | A | 0 | 0 | 0 | 0 | 100 | 40 | 100 | 23 | 36 | 47 | 58 | 55 | 57 | 42 |
| 1 | B | 0 | 0 | 0 | 0 | 100 | 15 | 100 | 16 | 9 | 38 | 51 | 46 | 54 | 37 |
| 2 | A | 0 | 0 | 0 | 0 | 57 | 12 | 0 | 0 | 10 | 36 | 14 | 14 | 25 | 21 |
| 2 | B | 0 | 0 | 0 | 0 | 18 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 24 | 19 |
| 3 | A | 1 | 0 | 0 | 0 | 100 | 88 | 100 | 12 | 76 | 66 | 71 | 65 | 72 | 76 |
| 3 | B | 1 | 0 | 0 | 0 | 100 | 89 | 100 | 11 | 75 | 67 | 72 | 65 | 72 | 75 |
| 4 | A | 1 | 0 | 0 | 0 | 100 | 78 | 100 | 4 | 78 | 62 | 69 | 60 | 58 | 63 |
| 4 | B | 1 | 0 | 0 | 0 | 100 | 95 | 100 | 15 | 57 | 69 | 73 | 68 | 54 | 78 |
| 5 | A | 0 | 0 | 1 | 0 | 100 | 90 | 100 | 0 | 47 | 67 | 72 | 62 | 47 | 71 |
| 5 | B | 0 | 0 | 1 | 0 | 100 | 0 | 100 | 0 | 0 | 32 | 46 | 37 | 68 | 22 |
| 6 | A | 0 | 0 | 1 | 0 | 100 | 59 | 100 | 16 | 35 | 55 | 63 | 59 | 53 | 73 |
| 6 | B | 0 | 0 | 1 | 0 | 100 | 18 | 100 | 6 | 14 | 39 | 52 | 44 | 53 | 22 |
| 7 | A | 1 | 0 | 1 | 0 | 100 | 34 | 100 | 0 | 25 | 45 | 56 | 46 | 70 | 59 |
| 7 | B | 1 | 0 | 1 | 0 | 100 | 100 | 100 | 67 | 60 | 71 | 75 | 86 | 64 | 79 |
| 8 | A | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 25 | 7 |
| 8 | B | 1 | 0 | 1 | 0 | 100 | 100 | 100 | 100 | 100 | 71 | 75 | 96 | 84 | 92 |
| 9 | A | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 14 | 8 |
| 9 | B | 0 | 1 | 0 | 0 | 100 | 100 | 100 | 0 | 100 | 71 | 75 | 65 | 78 | 90 |
| 10 | A | 0 | 1 | 0 | 0 | 77 | 18 | 0 | 0 | 22 | 39 | 15 | 16 | 15 | 23 |
| 10 | B | 0 | 1 | 0 | 0 | 98 | 79 | 0 | 0 | 63 | 63 | 33 | 33 | 21 | 55 |
| 11 | A | 1 | 1 | 0 | 0 | 100 | 70 | 100 | 77 | 68 | 59 | 66 | 80 | 71 | 93 |
| 11 | B | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 16 | 4 |
| 12 | A | 1 | 1 | 0 | 0 | 100 | 82 | 100 | 74 | 83 | 64 | 70 | 83 | 57 | 77 |
| 12 | B | 1 | 1 | 0 | 0 | 100 | 0 | 100 | 12 | 24 | 32 | 46 | 41 | 53 | 37 |
| 13 | A | 0 | 1 | 1 | 0 | 67 | 34 | 0 | 0 | 35 | 45 | 20 | 20 | 14 | 8 |
| 13 | B | 0 | 1 | 1 | 0 | 70 | 35 | 0 | 0 | 35 | 46 | 20 | 20 | 21 | 64 |
| 14 | A | 0 | 1 | 1 | 0 | 97 | 91 | 0 | 0 | 59 | 67 | 36 | 36 | 21 | 22 |
| 14 | B | 0 | 1 | 1 | 0 | 91 | 77 | 0 | 0 | 51 | 62 | 32 | 32 | 20 | 18 |
| 15 | A | 1 | 1 | 1 | 0 | 100 | 1 | 100 | 99 | 10 | 32 | 47 | 68 | 71 | 76 |
| 15 | B | 1 | 1 | 1 | 0 | 100 | 1 | 100 | 100 | 10 | 32 | 47 | 68 | 71 | 76 |
| 16 | A | 1 | 1 | 1 | 0 | 100 | 23 | 100 | 100 | 35 | 41 | 53 | 74 | 80 | 92 |
| 16 | B | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 14 | 4 |
| 17 | A | 0 | 0 | 0 | 1 | 100 | 19 | 100 | 37 | 18 | 39 | 52 | 54 | 66 | 69 |
| 17 | B | 0 | 0 | 0 | 1 | 100 | 0 | 100 | 36 | 17 | 32 | 46 | 48 | 65 | 46 |
| 18 | A | 0 | 0 | 0 | 1 | 100 | 11 | 100 | 40 | 17 | 36 | 49 | 52 | 55 | 63 |
| 18 | B | 0 | 0 | 0 | 1 | 100 | 7 | 100 | 37 | 9 | 34 | 48 | 50 | 56 | 66 |
| 19 | A | 1 | 0 | 0 | 1 | 100 | 74 | 100 | 7 | 65 | 61 | 67 | 60 | 55 | 53 |
| 19 | B | 1 | 0 | 0 | 1 | 100 | 72 | 100 | 7 | 65 | 60 | 67 | 59 | 55 | 49 |
| 20 | A | 1 | 0 | 0 | 1 | 100 | 92 | 100 | 8 | 72 | 68 | 72 | 65 | 57 | 41 |
| 20 | B | 1 | 0 | 0 | 1 | 100 | 88 | 100 | 4 | 53 | 66 | 71 | 63 | 56 | 71 |
| 21 | A | 0 | 0 | 1 | 1 | 100 | 47 | 100 | 40 | 45 | 50 | 60 | 63 | 58 | 80 |
| 21 | B | 0 | 0 | 1 | 1 | 100 | 9 | 100 | 21 | 10 | 35 | 49 | 46 | 59 | 18 |
| 22 | A | 0 | 0 | 1 | 1 | 100 | 100 | 100 | 89 | 83 | 71 | 75 | 92 | 48 | 60 |
| 22 | B | 0 | 0 | 1 | 1 | 100 | 8 | 100 | 0 | 15 | 35 | 49 | 39 | 53 | 42 |
| 23 | A | 1 | 0 | 1 | 1 | 5 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 15 | 3 |
| 23 | B | 1 | 0 | 1 | 1 | 91 | 79 | 0 | 0 | 72 | 63 | 33 | 33 | 22 | 39 |
| 24 | A | 1 | 0 | 1 | 1 | 100 | 66 | 100 | 4 | 63 | 58 | 65 | 57 | 57 | 44 |
| 24 | B | 1 | 0 | 1 | 1 | 100 | 94 | 100 | 22 | 79 | 69 | 73 | 70 | 54 | 73 |
| 25 | A | 0 | 1 | 0 | 1 | 100 | 25 | 100 | 21 | 26 | 42 | 53 | 50 | 69 | 43 |
| 25 | B | 0 | 1 | 0 | 1 | 100 | 74 | 100 | 54 | 65 | 61 | 67 | 74 | 56 | 73 |
| 26 | A | 0 | 1 | 0 | 1 | 100 | 6 | 100 | 3 | 9 | 34 | 48 | 40 | 60 | 39 |
| 26 | B | 0 | 1 | 0 | 1 | 100 | 87 | 100 | 35 | 54 | 66 | 71 | 72 | 46 | 69 |
| 27 | A | 1 | 1 | 0 | 1 | 100 | 97 | 100 | 52 | 97 | 70 | 74 | 80 | 67 | 80 |
| 27 | B | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 17 | 6 |
| 28 | A | 1 | 1 | 0 | 1 | 100 | 90 | 100 | 22 | 80 | 67 | 72 | 69 | 79 | 89 |
| 28 | B | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 12 | 5 |
| 29 | A | 0 | 1 | 1 | 1 | 100 | 58 | 100 | 24 | 44 | 55 | 63 | 61 | 66 | 47 |
| 29 | B | 0 | 1 | 1 | 1 | 100 | 63 | 100 | 24 | 38 | 57 | 64 | 62 | 54 | 67 |
| 30 | A | 0 | 1 | 1 | 1 | 100 | 57 | 100 | 29 | 49 | 54 | 63 | 62 | 63 | 58 |
| 30 | B | 0 | 1 | 1 | 1 | 100 | 46 | 100 | 24 | 39 | 50 | 60 | 57 | 63 | 56 |
| 31 | A | 1 | 1 | 1 | 1 | 100 | 2 | 100 | 4 | 3 | 32 | 47 | 39 | 63 | 44 |
| 31 | B | 1 | 1 | 1 | 1 | 100 | 4 | 100 | 4 | 4 | 33 | 47 | 39 | 53 | 46 |
| 32 | A | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 32 | 10 | 10 | 16 | 5 |
| 32 | B | 1 | 1 | 1 | 1 | 100 | 75 | 100 | 66 | 73 | 61 | 68 | 78 | 65 | 71 |
4.5.3.6 Heuristic model
fit.brm_exp3_causal_heuristic %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(term = str_remove(term, "b_"),
term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
# print_table(format = "latex")
print_table()| term | estimate | std.error | lower 95% HDI | upper 95% HDI |
|---|---|---|---|---|
| intercept | 49.79 | 1.46 | 47.39 | 52.23 |
| moving | 0.22 | 0.22 | 0.01 | 0.66 |
| speed | 2.06 | 0.85 | 0.64 | 3.43 |
| contact_e | 0.38 | 0.35 | 0.02 | 1.07 |
| e_speed_diff | 0.12 | 0.12 | 0.01 | 0.37 |
| e_direction_diff | 1.04 | 0.74 | 0.08 | 2.43 |
| total_speed_diff | 2.21 | 0.95 | 0.67 | 3.80 |
| total_direction_diff | 3.99 | 0.91 | 2.45 | 5.45 |
| transfer | 15.59 | 0.79 | 14.30 | 16.90 |
| e_moving | 0.18 | 0.17 | 0.01 | 0.51 |
| exclusive | 4.39 | 0.73 | 3.16 | 5.56 |
5 Appendix
5.0.0.1 Experiment 3
5.0.0.1.1 Does the question framing matter? (Responsibility vs. causation)
We have run a pilot eye-tracking experiment with 15 participants who saw the same clips as the ones in Experiment 3. Instead of asking participants to judge “To what extent were A and B responsible for E (not) going through the gate?”, we asked them to judge to what extent they agreed with the following statement: “Ball A/B caused ball E to go through the gate.” or “Ball A/B prevented ball E from going through the gate.” depending on the outcome.
Participants agreement judgments with the causal statements in this pilot eye-tracking experiment were strikingly similar to participants’ responsibility judgments in the online experiment.
df.exp3.eye_tracking_pilot = read_csv("../../data/experiment3_eye_tracking_pilot.csv") %>%
rename(clip = trial) %>%
mutate(ball = str_to_upper(ball)) %>%
group_by(clip, ball, outcome) %>%
summarize(causality_mean = mean(causation),
causality_low = smean.cl.boot(causation)[2],
causality_high = smean.cl.boot(causation)[3]) %>%
ungroup()
df.plot = df.exp3.eye_tracking_pilot %>%
# select(clip, ball, outcome,) %>%
left_join(df.exp3.causal.means %>%
select(clip,
ball,
responsibility_mean = rating_mean,
responsibility_low = rating_low,
responsibility_high = rating_high),
by = c("clip", "ball"))
ggplot(data = df.plot,
mapping = aes(x = causality_mean,
y = responsibility_mean)) +
geom_abline(intercept = 0,
slope = 1,
linetype = 2) +
geom_linerange(mapping = aes(ymin = responsibility_low,
ymax = responsibility_high),
alpha = 0.1) +
geom_linerange(mapping = aes(xmin = causality_low,
xmax = causality_high),
alpha = 0.1) +
geom_point(size = 2) +
geom_smooth(method = "lm",
color = "black") +
annotate(geom = "text",
x = c(0, 0),
y = c(100, 90),
label = c(str_c("r = ",
cor(df.plot$causality_mean,df.plot$responsibility_mean) %>%
round(2)),
str_c("RMSE = ", rmse(df.plot$causality_mean,df.plot$responsibility_mean) %>%
round(2))),
hjust = 0,
size = 8) +
scale_x_continuous(breaks = seq(0, 100, 25),
limits = c(0, 100)) +
scale_y_continuous(breaks = seq(0, 100, 25),
limits = c(0, 100)) +
labs(x = "causal judgment",
y = "responsibility judgment") +
theme(text = element_text(size = 24))
# ggsave("../../figures/plots/aux_question_framing_comparison.pdf",
# width = 6,
# height = 6)
5.0.0.1.2 3D scatter plot of participant clusters
plot_ly(x = df.cluster_whs$whether,
y = df.cluster_whs$how,
z = df.cluster_whs$sufficient,
type = "scatter3d",
mode = "markers",
color = as.factor(df.cluster_whs$cluster),
colors = c("#e41a1c", "#377eb8")) %>%
layout(showlegend = FALSE,
title = "Participant clusters",
scene = list(
xaxis = list(title = "whether"),
yaxis = list(title = "how"),
zaxis = list(title = "sufficient")))5.0.0.1.3 Model fits for individual participants
df.plot = fit.brm_exp3_causal_whs %>%
fitted() %>%
as_tibble() %>%
select(prediction = Estimate) %>%
bind_cols(df.exp3.causal.long) %>%
relocate(prediction, .after = last_col())
df.text = df.plot %>%
group_by(participant) %>%
summarize(r = round(cor(prediction, rating), 2)) %>%
ungroup() %>%
mutate(r = str_c("r = ", r),
prediction = 1,
rating = 110)`summarise()` ungrouping output (override with `.groups` argument)ggplot(data = df.plot,
mapping = aes(x = prediction,
y = rating)) +
geom_abline(intercept = 0,
slope = 1,
color = "blue",
alpha = 0.5) +
geom_point(alpha = 0.3,
size = 1) +
geom_text(data = df.text,
mapping = aes(label = r),
hjust = 0,
color = "red",
size = 4) +
facet_wrap(facets = vars(participant),
ncol = 7) +
labs(x = "model prediction",
y = "causal responsibility rating") +
coord_cartesian(xlim = c(0, 100),
ylim = c(0, 100),
expand = F,
clip = "off") +
scale_x_continuous(breaks = seq(0, 100, 25),
limits = c(0, 100)) +
scale_y_continuous(breaks = seq(0, 100, 25),
expand = expansion(mult = c(0, 0))) +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_blank(),
text = element_text(size = 16),
panel.spacing.x = unit(0.75, "cm"),
panel.spacing.y = unit(1, "cm"),
plot.margin = margin(t = 0.5, r = 0.5, unit = "cm"))
# ggsave("../../figures/plots/exp3_individual_scatter_plots.pdf",
# width = 12,
# height = 8)
5.0.0.1.4 Individual participant model fits
df.plot = df.exp3.causal.individual_fit %>%
select(participant, contains("r_"), contains("rmse_")) %>%
pivot_longer(cols = -participant,
names_to = c("measure", "model"),
names_sep = "_",
values_to = "fit")
ggplot(data = df.plot %>%
filter(measure == "r"),
mapping = aes(x = fit,
fill = model)) +
geom_density(alpha = 0.5) +
labs(x = "correlation value") +
scale_x_continuous(breaks = seq(0, 1, 0.2),
limits = c(0, 1)) +
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
scale_fill_brewer(palette = "Set1") +
# ggplot2::theme_classic() +
theme(legend.position = c(0.3, 0.95),
legend.direction = "horizontal",
text = element_text(size = 20))
# ggsave("../../figures/plots/exp3_clusters_densities.pdf",
# width = 8,
# height = 6)
5.0.0.1.5 Clusters of participants
df.plot = df.exp3.causal.long %>%
left_join(df.cluster_whs %>%
select(participant, cluster) %>%
mutate(participant = as.numeric(participant),
cluster = as.factor(cluster)) %>%
group_by(cluster) %>%
mutate(label = factor(cluster, levels = 1:2, labels = str_c("n = ", group_size(.)))) %>%
ungroup(),
by = "participant")
ggplot(data = df.plot,
mapping = aes(x = ball,
y = rating,
group = label,
fill = label)) +
geom_point(shape = 21,
alpha = 0.2,
position = position_jitterdodge(jitter.width = 0.1,
jitter.height = 0,
dodge.width = 0.75)) +
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
mapping = aes(color = label),
size = 1,
position = position_dodge(width = 0.75)) +
stat_summary(fun = "mean",
geom = "point",
size = 4,
shape = 21,
position = position_dodge(width = 0.75)) +
facet_wrap(facets = vars(clip), ncol = 8) +
labs(y = "causal responsibility rating") +
scale_fill_brewer(palette = "Set1") +
scale_color_brewer(palette = "Set1") +
theme(text = element_text(size = 24),
legend.position = "bottom")
# ggsave("../../figures/plots/exp3_clusters_points.pdf",
# width = 12,
# height = 8)
5.0.0.1.6 Ternary plots for individual regression results
Registered S3 methods overwritten by 'ggtern':
method from
grid.draw.ggplot ggplot2
plot.ggplot ggplot2
print.ggplot ggplot2--
Remember to cite, run citation(package = 'ggtern') for further info.
--
Attaching package: 'ggtern'The following objects are masked from 'package:gridExtra':
arrangeGrob, grid.arrangeThe following objects are masked from 'package:ggplot2':
aes, annotate, ggplot, ggplot_build, ggplot_gtable, ggplotGrob,
ggsave, layer_data, theme_bw, theme_classic, theme_dark,
theme_gray, theme_light, theme_linedraw, theme_minimal, theme_voiddf.plot = df.exp3.causal.individual_fit %>%
select(participant, fit_whs) %>%
mutate(estimates = map(fit_whs, tidy)) %>%
select(participant, estimates) %>%
unnest(estimates) %>%
filter(str_detect(term, "b_"),
!str_detect(term, "Intercept")) %>%
mutate(term = str_remove(term, "b_")) %>%
select(participant, term, estimate) %>%
pivot_wider(names_from = term,
values_from = estimate) %>%
mutate_at(vars(whether, how, sufficient),
.funs = list(norm = ~ . / (how + whether + sufficient))) %>%
mutate(color = 0,
color = ifelse(how_norm == max(how_norm), 1, color),
color = ifelse(whether_norm == max(whether_norm), 2, color),
color = ifelse(sufficient_norm == max(sufficient_norm), 3, color),
color = factor(color))
ggplot(data = df.plot,
mapping = aes(x = how,
y = sufficient,
z = whether,
color = color)) +
geom_point(alpha = 0.7,
size = 2,
show.legend = F) +
scale_color_manual(values = c("black", "red", "green", "blue")) +
coord_tern() +
theme_showarrows() +
theme(text = element_text(size = 20),
tern.axis.title.T = element_blank(),
tern.axis.title.L = element_blank(),
tern.axis.title.R = element_blank())Coordinate system already present. Adding new coordinate system, which will replace the existing one.# ggsave(str_c("../../figures/plots/exp3_individual_regression_ternary_plot_scaled.pdf"),
# width = 5,
# height = 5)
6 Session info
R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] ggtern_3.3.0 forcats_0.5.0 stringr_1.4.0 dplyr_1.0.0
[5] purrr_0.3.4 readr_1.3.1 tidyr_1.1.0 tibble_3.0.1
[9] tidyverse_1.3.0 patchwork_1.0.0 egg_0.4.5 gridExtra_2.3
[13] plotly_4.9.2.1 png_0.1-7 tidybayes_2.0.3 Hmisc_4.4-0
[17] ggplot2_3.3.1 Formula_1.2-3 survival_3.1-12 lattice_0.20-41
[21] janitor_2.0.1 jsonlite_1.6.1 xtable_1.8-4 corrr_0.4.2
[25] brms_2.13.0 Rcpp_1.0.4.6 broom_0.5.6 lme4_1.1-23
[29] Matrix_1.2-18 kableExtra_1.1.0 knitr_1.28
loaded via a namespace (and not attached):
[1] proto_1.0.0 tidyselect_1.1.0 htmlwidgets_1.5.1
[4] munsell_0.5.0 codetools_0.2-16 statmod_1.4.34
[7] DT_0.13 miniUI_0.1.1.1 withr_2.2.0
[10] Brobdingnag_1.2-6 colorspace_1.4-1 highr_0.8
[13] rstudioapi_0.11 stats4_3.6.2 robustbase_0.93-6
[16] bayesm_3.1-4 bayesplot_1.7.2 labeling_0.3
[19] emmeans_1.4.7 rstan_2.19.3 farver_2.0.3
[22] bridgesampling_1.0-0 coda_0.19-3 vctrs_0.3.1
[25] generics_0.0.2 TH.data_1.0-10 xfun_0.14
[28] R6_2.4.1 markdown_1.1 HDInterval_0.2.2
[31] assertthat_0.2.1 promises_1.1.0 scales_1.1.1
[34] multcomp_1.4-13 nnet_7.3-14 gtable_0.3.0
[37] processx_3.4.2 sandwich_2.5-1 rlang_0.4.6
[40] splines_3.6.2 lazyeval_0.2.2 acepack_1.4.1
[43] checkmate_2.0.0 inline_0.3.15 yaml_2.2.1
[46] reshape2_1.4.4 abind_1.4-5 modelr_0.1.8
[49] threejs_0.3.3 crosstalk_1.1.0.1 backports_1.1.7
[52] httpuv_1.5.3.1 rsconnect_0.8.16 tensorA_0.36.1
[55] tools_3.6.2 bookdown_0.19 ellipsis_0.3.1
[58] RColorBrewer_1.1-2 ggridges_0.5.2 latex2exp_0.4.0
[61] plyr_1.8.6 base64enc_0.1-3 ps_1.3.3
[64] prettyunits_1.1.1 rpart_4.1-15 zoo_1.8-8
[67] haven_2.3.1 cluster_2.1.0 fs_1.4.1
[70] magrittr_1.5 data.table_1.12.8 colourpicker_1.0
[73] reprex_0.3.0 mvtnorm_1.1-0 matrixStats_0.56.0
[76] hms_0.5.3 shinyjs_1.1 mime_0.9
[79] evaluate_0.14 arrayhelpers_1.1-0 shinystan_2.5.0
[82] jpeg_0.1-8.1 readxl_1.3.1 rstantools_2.0.0
[85] compiler_3.6.2 crayon_1.3.4 minqa_1.2.4
[88] StanHeaders_2.21.0-3 htmltools_0.4.0 mgcv_1.8-31
[91] later_1.1.0.1 lubridate_1.7.8 DBI_1.1.0
[94] dbplyr_1.4.4 MASS_7.3-51.6 boot_1.3-25
[97] compositions_1.40-5 cli_2.0.2 parallel_3.6.2
[100] igraph_1.2.5 pkgconfig_2.0.3 foreign_0.8-76
[103] xml2_1.3.2 svUnit_1.0.3 dygraphs_1.1.1.6
[106] webshot_0.5.2 estimability_1.3 rvest_0.3.5
[109] snakecase_0.11.0 callr_3.4.3 digest_0.6.25
[112] rmarkdown_2.2 cellranger_1.1.0 htmlTable_1.13.3
[115] shiny_1.4.0.2 gtools_3.8.2 nloptr_1.2.2.1
[118] lifecycle_0.2.0 nlme_3.1-148 viridisLite_0.3.0
[121] fansi_0.4.1 pillar_1.4.4 loo_2.2.0
[124] fastmap_1.0.1 httr_1.4.1 DEoptimR_1.0-8
[127] pkgbuild_1.0.8 glue_1.4.1 xts_0.12-0
[130] shinythemes_1.1.2 stringi_1.4.6 blob_1.2.1
[133] latticeExtra_0.6-29